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Abstract 

The focus of the current study is on handling the dependence among multiple regression coefficients representing the treatment 
effects when meta-analyzing data from single-case experimental studies. We compare the results when applying three different 
multilevel meta-analytic models (1.¢., a univariate multilevel model avoiding the dependence, a multivariate multilevel model 
ignoring covariance at higher levels, and a multivariate multilevel model modeling the existing covariance) to deal with the 
dependent effect sizes. The results indicate better estimates of the overall treatment effects and variance components when a 
multivariate multilevel model is applied, independent of modeling or ignoring the existing covariance. These findings confirm 
the robustness of multilevel modeling to misspecifying the existing covariance at the case and study level in terms of estimating 
the overall treatment effects and variance components. The results also show that the overall treatment effect estimates are 
unbiased regardless of the underlying model, but the between-case and between-study variance components are biased in certain 
conditions. In addition, the between-study variance estimates are particularly biased when the number of studies is smaller than 
40 (i.e., 10 or 20) and the true value of the between-case variance is relatively large (i.e., 8). The observed bias is larger for the 
between-case variance estimates compared to the between-study variance estimates when the true between-case variance is 
relatively small (1.e., 0.5). 


Keywords Single-case experimental design - Multilevel meta-analysis - Multivariate multilevel model - Robust variance estimator 


Introduction 


Single-case experimental designs (SCEDs) have been fre- 
quently used in many different fields such as psychology, 
education, social sciences, and medical sciences to establish 
evidence-based practices. In this kind of design, the outcome 
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of interest is measured repeatedly over time prior to and dur- 
ing or after the treatment. Afterwards, scores under these con- 
ditions (pre and post treatment) can be compared to evaluate 
treatment effectiveness (Horner et al., 2005; Kratochwill et al., 
2010). An inherent SCED characteristic is that each identity 
(e.g., an individual or a group of individuals) serves as its own 
control. There are a variety of SCED types including with- 
drawal or reversal designs, alternating treatment designs, 
changing criterion designs, and multiple baseline designs 
(MBDs). In the current study, we focus on MBDs in which 
the introduction of the treatment is staggered across different 
time points for multiple cases, outcomes, or settings. The rea- 
son for this choice was that MBDs are not only the most 
frequently used SCED design, but they are also internally 
and externally more valid compared to other SCED design 
types (Kratochwill et al., 2010; Onghena & Edgington, 2005). 

The data from several MBDs can be synthesized by apply- 
ing multilevel meta-analytic modeling proposed by Van den 
Noortgate and Onghena (2003a, 2003b, 2008). A three-level 
model can be used to model the hierarchical structure of the 
data: the measurement occasions (i.e., the first-level units) are 


a Springer 


nested within cases (1.e., the second-level units), which in turn 
are nested within studies (i.e., the third-level units). Multilevel 
modeling not only can yield improved estimates of the treat- 
ment effects for individual cases, but also an average treatment 
effect estimate across cases for each study and an overall treat- 
ment effect estimate across all studies. The variance at the first 
level, that is the variance between measurement occasions 
within cases, can be modeled as follows: 


with ej,.~N (0, a) 
(1) 


In this equation, i is an indicator for the measurement oc- 
casion (i= 0, 1, ..., ),/ for the case (j = 0, 1, ..., J), and & for 
the study (k=0, 1, ..., K). Therefore, Y;, is the outcome on the 
i” measurement occasion from the j” case within study k. The 

dummy variable D,,; indicates the corresponding phase for 
each measurement and equals 0 if the measurement belongs 
to the baseline phase and | if it belongs to the treatment phase. 
Tix is a time indicator (e.g., days, the session number). If we 
center the time indicator such that it is equal to 0 at the start 
point of the treatment phase, /3;, expresses the expected base- 
line level at the start point of the treatment phase (where 7), = 
0) and (35;, the immediate treatment effect on the outcomes. 
The coefficient (31; refers to the linear time trend in the out- 
come during the baseline phase, and (33; indicates the treat- 
ment effect on the time trend. The sampling error, e;x, is as- 
sumed normally, identically and independently distributed 
with a zero mean and a sampling variance of o2. The indices 
j and k in the regression coefficients mean that these coeffi- 
cients are case- and study-specific, respectively. The variation 
over cases (second level) is modeled using additional equa- 
tions. In a simple model, the case-specific regression coeffi- 
cients are assumed to deviate randomly from the study mean 
coefficients. The following equations show that the baseline 
level, time effect in the baseline phase, the immediate treat- 
ment effect, and the treatment effect on the time trend, respec- 
tively, vary randomly around the study-specific means: 


Vie = Boje + Pry T ie + BryeDie + Ban T nD + eit 


Bojx = Pook + Uojk Udjk 
Bre = Prox + Uryx Ul; 
: jk 
as z with | ~'* |~MVN(0,»,) 
Boje = 920K + Urjx U2; 
Bsjn = 930% + U3yx U3; 
2 
Tx 4 
—_ Oulud On 
Ly = 2 
Ou2ud0 = AFu2ul O12 
Oy3u0 = =Fu3ul = =Fu3u2 os 


(2) 


Each case-specific coefficient 3 jx (1-€., Boje, Pijks Boje. OF 
(jx) equals a study-specific mean 4 ox (1.€., Poors Prox F20% 
or 430,) plus a random case-specific deviation from this mean. 
The covariance matrix , has on its diagonal variances 
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between cases within the same study andand the off- 
diagonal elements the covariances between cases within the 
same study. For instance, 07, and 02, are the variances be- 
tween cases from the same study of the immediate treatment 
effect and the treatment effect on the time trend; 0,,3,,2 repre- 
sents the covariance between the immediate treatment effect 
and the treatment effect on the time trend. 

The variation of the study-specific means, @’s, are modeled 
at the third level as the overall regression coefficients, ’s, plus 
a random study-specific deviation from these overall means 
over all studies. 


Pook = Yooo + Voor Vook 
0 = +r Vv = Vv 
10k  "Yioo F V10k with | “0K | ~MVN(0, 5,) 
920k = Y200 + V20K V20k 
30x = Y300 + V30K V30k 
2 
oo > 
xy, = Ov1v0 1 
a o 
Oy2v0  Ay2vl v2 i 
130 = F3vl =Fy3v2 3 
(3) 


Residuals at case and study levels are assumed to be multi- 
variate normally distributed. Researchers are mainly interested 
in the overall treatment effects (i.e., the overall immediate treat- 
ment effect, 299, and the overall treatment effect on time trend, 
7300) and the (co)variation in these effects over cases and stud- 
ies. These (and other) model parameters can be estimated by 
applying maximum likelihood procedures commonly used in 
multilevel modeling. In addition, one can also obtain improved 
estimates of the case- and study-specific treatment effects (the 
(’s and 6’s) by using empirical Bayes techniques (Raudenbush 
& Bryk, 2002). Multilevel modeling has been applied in the 
context of SCED data in previous studies, and its performance 
has been assessed in several simulation studies (e.g., Moeyaert, 
Ugille, Ferron, Beretvas, & Van den Noortgate, 2013; Owens & 
Ferron, 2012; Ugille, Moeyaert, Beretvas, Ferron, & Van den 
Noortgate, 2012, 2014). In general, the findings have indicated 
appropriate performance in terms of statistical properties such 
as relative bias, MSE of this method, both for synthesizing the 
standardized and unstandardized SCED data. 

The model can be modified in various ways, according to 
the data characteristics and the research interests. For instance, 
predictor variables can be added at each level to study their 
effect on the outcome, or to study how they moderate the 
treatment effects. In addition, it is possible to use alternative 
specifications of the structure of variances and covariances at 
each level. 

The model can also be adapted to combine effect sizes, 
rather than raw data from multiple SCED studies. If we obtain 
the estimates for two effect sizes of interest (i.e., the immediate 
treatment effect and the treatment effect on the time trend) for 
each case separately by fitting a simple linear model (Eq. 1), 
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we can then synthesize the case-specific ordinary least squares 
(OLS) estimates (i.¢., bz; and b3;,) by using either a separate 
univariate three-level meta-analysis or a multivariate three- 
level meta-analysis (Van den Noortgate & Onghena, 2003b, 
2008). For instance, by; is equal to 2;,, the immediate treat- 
ment effect for the case 7 from the study k, plus a random 
deviation that is due to sampling variance. At the second and 
third levels, the equations for immediate effects are the same 
as for the raw data model (see Eqs. 2 and 3). We can rewrite 
this three-level model in a single equation (Eq. 4): 


bojk = 200 + Urjk + V20K + 12k 


rajp~N (0. oy) > urjp~N (0. 7.) ; v2jnh~N (0. o,) (4) 

in which this estimate of the immediate effect for case j 
from study k is equal to the overall immediate treatment effect, 
200. plus a case-specific random deviation, u2;;,, a study- 
specific random deviation, v29,, and a residual deviation due 
to sampling variance, r;;. This model provides the overall 
treatment effect, 290, as well as the corresponding between- 
case and between-study variance components, oj, and oy, 
respectively. From the initial OLS regression analysis, we al- 
ready have an estimate of the sampling variance, on , So in the 
meta-analysis it is assumed to be known. We can use the same 
model for synthesizing the estimates of the treatment effect on 
the time trend, b3,, by using Eq. 5: 


b3jk = Y300 + U3ik + V30K + 13k 


ryx~N (0, oy.) 9 u3;k~N (0. o,) 5 V3;h~N (0. o,) (5) 

Alternatively, another way to analyze the multiple effect 
sizes, taking into account the dependence among these effect 
sizes, is the use of the multivariate mixed model proposed by 
Kalaian and Raudenbush (1996), Eq. 6. This model uses the 
estimated sampling covariance matrix of multivariate effect 
sizes to obtain the optimal weights for the effect sizes and to 
estimate model parameters and their corresponding standard 
errors rather than using sampling variance estimates of the 
effect sizes, which is the case in ordinary meta-analysis. One 
of the most important advantages of this approach is that all 
available information is used in a single analysis in order to 
estimate any treatment effects (Van den Noortgate, Lopez- 
Lopez, Marin-Martinez, & Sanchez-Meca, 2013). 


bik = Y200 + Uajk + V20K + Tr; 
bse = Y300 + Usjk + V30K + 135k 


(2%) mewn (0,2°(0)), (2" )~aeviv(a, 335), ( 72" )~a0viv(0, 


T3 ik U3;k 3jk 
(6) 


The sampling covariance matrix, ’(b), the between-case 
and between-study covariance matrices 41,, and 4, are elements 


of R°*?. Combining these case-specific effect sizes using a 
bivariate three-level meta-analytic model, Eq. 6, results in the 
estimates of two fixed effects (7209 and 7399) and six unique 
(co)variance components (the elements of 4, and },). 

The main advantage of the multivariate model over sepa- 
rate univariate three-level models is its ability to provide more 
information about the potential covariances between effect 
sizes. This strategy explicitly models the correlations among 
the effect size estimates and needs information regarding the 
covariance structure of the effect size estimation errors. 
Getting knowledge about this structure requires the informa- 
tion about the dependence structure of the raw data in each 
study (Hedges, Tipton, & Johnson, 2010). If the raw data are 
available, an estimate of the sampling covariance between 
effect sizes can indeed be obtained from the OLS regression. 
A major drawback of this multivariate approach is that when 
the raw data are not available, this approach is often not ap- 
plicable because it is often possible to reconstruct the sam- 
pling variances (the squared standard errors, that are typically 
reported or can be reconstructed based on the results of statis- 
tical significance tests), but not the sampling covariances (Van 
den Noortgate et al., 2013). 

To understand the sampling covariance between both ef- 
fects (1.e., immediate treatment effect and treatment effect on 
the time trend) at the lowest level, consider the situation where 
there is a positive time trend in the baseline and a positive 
immediate treatment effect and a positive effect on the time 
trend. If by coincidence the first measurement in the treatment 
phase happens to be unexpectedly large, the immediate effect 
is likely to be overestimated, while the effect on the trend will 
be underestimated. In general, a negative correlation between 
both coefficients can be expected. But also at the level of 
participants or studies, both kinds of effect might be depen- 
dent: for instance, both coefficients can be negatively corre- 
lated at the case level, if, for example, some cases show im- 
mediately a full effect of the treatment, while for other cases, 
the effect is building up gradually. In this situation, a non-zero 
covariance between residuals at the case and study level can 
be expected. Still, in most previous studies of multilevel 
modeling to synthesize the standardized and unstandardized 
SCED data (e.g., Moeyaert et al., 2013; Ugille et al., 2012, 
2014), it was assumed that each random effect is independent- 
ly, identically, and normally distributed with the mean of zero 
and a diagonal covariance matrix at the second and third level. 
Moeyaert, Ugille, Ferron, Beretvas, and Van den Noortgate 
(2016) investigated the impact of covariance matrix 
misspecification on the quality of the estimates of fixed and 
random effects in the situation where the residuals covary at 
the case and study level. They modeled the condition in which 
there were two regression coefficients (1.e., the baseline level 
and the treatment level) in the multilevel modeling. They 
found that misspecifying the existing covariance has no sig- 
nificant impact on the estimates of treatment effect. Moreover, 
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the results indicated that modeling or ignoring the covariance 
leads to relatively unbiased between-case variance estimates, 
whereas the between-study variance estimates are biased in 
both situations. 

The aim of current study is to compare different approaches 
to deal with dependent effect sizes in SCED research. We 
focus on the situation where you have two effect sizes per 
case: one referring to the immediate treatment effect and one 
referring to the effect on a time trend. To accommodate these 
dependent effect sizes, we use three different approaches: a 
traditional univariate random effects approach (for each type 
of coefficients separately) avoiding the existing dependence, a 
multivariate random effects meta-analytic approach proposed 
by Kalaian and Raudenbush (1996) ignoring the existing co- 
variance at higher levels, and a multivariate multilevel meta- 
analysis modeling the covariance at higher levels. In the last 
two approaches, we also study the performance of the sand- 
wich estimator (Huber, 1967; White, 1982), also known as the 
robust covariance matrix estimator, the empirical covariance 
matrix estimator, or robust variance estimator (RVE; Hedges 
et al., 2010). This procedure is assumed to yield asymptotical- 
ly consistent covariance matrix estimates, even if the underly- 
ing models are incorrect or misspecified. We would like to 
understand whether treating the existing dependence in differ- 
ent ways has any significant impact on the quality of the over- 
all treatment effects estimates and the between-case and 
between-study (co)variance estimates. A difference with the 
study of Moeyaert et al. (2016), is that we focus on the meta- 
analysis of effect sizes rather than on the meta-analysis of 
SCED raw data, and we compare the results of a multivariate 
approach with those of univariate analyses. A difference with 
the studies of Ugille et al. (2012, 2014), is that we explicitly 
assume that effect sizes are correlated at higher levels and that 
we do not only use univariate models to analyze the data. 


Method 


A simulation study was conducted in order to handle the 
dependence among multiple regression coefficients in the 
multilevel meta-analysis of SCED studies when freely es- 
timating covariances between residuals at higher levels 
(i.e., case-level and study-level). Several models were 
fitted to the generated data to evaluate their performance 
in parameter recovery. 

Data generation and analyses were done using SAS 9.4 
(Fan, Felsovalyi, Sivo, & Keenan, 2002; SAS Institute Inc., 
2017). 


Data generation 


SCED raw data were generated using the three-level model 
from Eqs. | to 3 for MBDs across participants. The 
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coefficients of covariates in the baseline phase (1. e. , Yooo 
and ‘y;99) were set to 0, whereas the coefficients of covariates 
related to the treatment (i.¢., 7209 and 7300) were equal to 2 and 
0.2, respectively. We manipulated the between-cases variabil- 
ity and the between-studies variability. In addition, we varied 
the size of simulated datasets by varying the number of mea- 
surement occasions J (J = 10, 20, or 40), the number of cases 
within studies J (J = 3, 4, or 7), and the number of studies 
meta-analyzed K (K = 10, 20, or 40). The values for J were 
chosen based on the findings of several meta-analyses and 
systematic reviews of SCEDs (Ferron, Farmer, & Owens, 
2010; Shadish & Sullivan, 2011). The values for J were se- 
lected based on recommendation of several studies on SCEDs 
(Farmer, Owens, Ferron, & Allsopp, 2010; Ferron et al., 2010; 
Moeyaert et al., 2013; Shadish & Sullivan, 2011) to find out 
how small the number of cases per study could be to return 
unbiased estimates in combination with other conditions. For 
K, we chose the values based on a systematic review of SCED 
meta-analyses (Jamshidi et al., 2020). 

Because the data were generated for MBDs, the timing of 
treatment introduction was staggered across subjects. Start 
points of the treatment for each case (Table 1) were chosen 
in line with Ugille et al. (2012) and Moeyaert, Ugille, Ferron, 
Beretvas, and Van den Noortgate (2014). 

Previous studies have shown that the between-study and 
between-case variances are significant factors affecting the 
estimates of treatment effects and the variances of these 
estimates (Moeyaert et al., 2013, 2014). The level-1 vari- 
ance of the raw data was set to | and assumed to be homo- 
geneous across phases. For the between-case and between- 
study covariance matrices (i.e., 4, and 4, the covariance 
matrices at the second and third level), we chose relatively 
small or relatively large amounts of variance, and therefore 
the variances (i.e., the diagonal elements in the matrices) 
were chosen to be either (8, 0.08, 8, 0.08) or (0.5, 0.05, 0.5, 
0.05). Because we were interested in the treatment effect 
estimates and not in the baseline levels, we solely assumed 
a non-zero covariance between the immediate treatment 
effect and treatment effect on the time trend at both levels 
(1.€., 0,243 and 02,3). We considered either no correlation 
(p = 0) or a large negative correlation (p = — .7, 


Table 1 Starting occasion of treatment introduction 
IJ 


3 (case | to case 3) 4 (case 1 to case 4) 7 (case | to case7) 


10 4,6,8 4,5,7,8 4, 5,5, 6, 7, 7,8 
20 7, 11,15 7, 10, 12, 15 79,9, UW, 13313,:15 
40 11, 21,31 11, 18, 24, 31 11,15, 15, 21,27, 27,31 


Notes: I = number of measurement occasions; J = number of cases 
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corresponding to a covariance of — 0.56 or — 0.11, depend- 
ing on the chosen population variances) between these 
treatment effects. This rather large correlation was chosen 
to make a possible effect of this factor clearly visible. The 
value chosen was negative in order to simulate the situa- 
tion described in the Introduction, where one may either 
show a larger immediate effect or a gradually increasing 
effect. Therefore, in &,, and §, all the off-diagonal ele- 
ments of the covariance matrices were set to 0, except for 
0243 aNd O,,2,3 for which the values of 0, — 0.56, or — 0.11 
were used. Crossing the values of the diagonal and off- 
diagonal elements led to four combinations for the 
between-case and four for the between-study covariance 
matrix (i.e., relatively large variance with no covariance, 
relatively large variance with covariance, relatively small 
variance with no covariance, and relatively small variance 
with covariance). 

In total, we evaluated 3 x 3x 3x 4x 4=432 conditions and 
for each combination, we simulated 2000 data sets, resulting 
in 864,000 data sets to be analyzed. 


Data analysis 


We fitted an OLS model to the generated raw data of each 
individual case, using Eq. 1, to estimate the regression coeffi- 
cients representing the treatment effect sizes (i.e., the immedi- 
ate effect and the effect on time trend, b;, and b3;,) and the 
residual’s standard deviation (Gj ). These case-specific inter- 
vention effect size estimates were then standardized by divid- 
ing them by the estimated residual standard deviation (Van 
den Noortgate & Onghena, 2003b). 


b; 
‘eas (7) 


For group-comparison designs, Hedges (1981) proposed to 
multiply observed standardized mean differences with an ap- 
proximate bias correction factor, in order to reduce the small 
sample bias from standardized mean differences: 


3 


ee (8) 


c(m) = 1 
with m equal to the degrees of freedom, which is the number 
of participants minus two. 

In a previous simulation study, we found that applying 
Hedges' (1981) bias correction to the standardized effect 
sizes leads to better estimates of overall treatment effects 
and variance components, especially for small / (Jamshidi 
et al., 2019). Therefore, in the current simulation, we syn- 
thesized the bias-corrected standardized regression coeffi- 
cients. For SCED data, the degrees of freedom, m, is equal 
to the number of measurement occasions, J, minus the 
number of predictors, p, in the regression model (Eq. 1) 


minus one (i.e., m = J — p - 1). According to Eq. 1, the 
degrees of freedom equals J - 4. The bias-corrected stan- 
dardized effect size is: 


(5) = 5 (1-5) (9) 


We also applied the bias correction factor to the standard 
error estimates and therefore to the sampling variance corre- 
sponding with each coefficient: 


Ae 2 5 3 2 
(#1) =23 (1-355) (10) 
We analyzed the bias-corrected standardized regression co- 


efficients, i.e., (04) and (sx) with three different models 


including a separate univariate three-level meta-analytic mod- 
el for each type of effect, a multivariate three-level meta-ana- 
lytic model ignoring the covariance structure at the case and 
study level, and a multivariate three-level meta-analytic model 
including the covariance structure at the case and study level. 
We used the last two models without and with RVE approach, 
so in total, five approaches were evaluated. In the following, 
“Model 1” refers to the two separate univariate three-level 
meta-analytic models for the immediate effect and the effect 
on time trend, using Eqs. 4 and 5. By analyzing the two types 
of effects separately, we avoid the issue of modeling depen- 
dence among these effect sizes in our analysis. 

In the second approach, we applied a multivariate three- 
level meta-analytic model (later referred to as “Model 2”), 
proposed by Kalaian and Raudenbush (1996), Eq. 6, but ig- 
noring the covariance structure at the case and study level. For 
doing so, first the Cholesky transformation (Kalaian & 
Raudenbush, 1996) was used to orthogonalize the regression 
coefficients for each case at the first level. Doing so, we could 
estimate the model parameters, using the sampling variance 
and covariance estimates obtained from the OLS regression 
analysis. A bivariate multilevel meta-analysis on these regres- 
sion coefficients was then used, estimating the overall effect 
on level and on the time trend, as well as their variances at the 
case and study level. In a third approach (later referred to as 
“Model 2 with RVE”), we used the sandwich estimator for the 
standard errors of the overall treatment effect estimates. The 
possible advantage of using RVE is that it can return asymp- 
totically consistent covariance matrix estimates, even if the 
underlying models are incorrect or misspecified (Huber, 
1967; White, 1982). We were interested in evaluating the 
functionality of adding this estimator to our model. RVE only 
adapts the standard errors, not the estimates of the fixed pa- 
rameters and variance components. Therefore, in the results 
section we only provide separate results for the second and 
third approach for the relative bias of standard errors 
estimates. 
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In the fourth approach, we used again a three-level multi- 
variate model, but this time including covariances at the case 
and study level (“Model 3”). In the fifth and final approach 
(referred to as “Model 3 with RVE”), we used the same model 
but applied the sandwich estimator to the standard errors esti- 
mates of the overall treatment effects. 

For all approaches, we made use of the restricted maximum 
likelihood procedure implement in PROC MIXED, using the 
Kenward-Roger approach to approximate the degrees of free- 
dom for the ¢-tests of the fixed effects. We selected this method 
for estimating degrees of freedom, because it was found to 
perform better in synthesizing multiple-baseline studies 
through multilevel modeling with a complex error structure 
(Ferron et al., 2010). 

To evaluate the performance of the five approaches for 
estimating and testing the overall treatment effects parameters 
(1.€., Y209 and ~y399), we used various criteria. The amount of 
absolute bias was approximated by subtracting the true value 
from the mean estimated value. We calculated the mean 
squared error (MSE) in order to check the accuracy of the 
estimates. Moreover, we looked at the estimates of the stan- 
dard error (SE) of the overall treatment effects’ estimates, that 
can be used to construct confidence intervals. The standard 
deviation of the treatment effects estimates can be assumed as 
an appropriate approximation of the true standard errors. The 
bias of the standard errors was therefore evaluated by the 
difference between the median standard error estimates and 
the standard deviation of the estimates of the treatment effects. 
Because the distributions of the standard errors tend to be 
skewed, we used the median rather than the mean of the stan- 
dard errors. The coverage proportions (CPs) of the 95% con- 
fidence intervals of the overall treatment effects estimates 
were constructed to evaluate the treatment effects estimates 
and the estimated standard errors. Given that for each condi- 
tion 2000 datasets are generated, the standard error of the 
estimate of the coverage proportion is approximately .005; 
therefore, a fairly good coverage proportion would be expect- 
ed to have a value between .94 and .96. An appropriate cov- 
erage proportion of the intervals would be a good indicator of 
a lack of bias in both the parameter estimates and the corre- 
sponding standard errors. 

We also looked at the performance of three-level models in 
estimating the variance components estimates ignoring and 
modeling the covariance at the case and study level. This 
way, we could evaluate the impact of the analysis approach 
on the between-case variance and between-study variance of 
the treatment effects in terms of accuracy of estimates (1.e., 
bias) and the precision of the estimates (i.e., MSE). 

Afterwards, we calculated the relative bias of each param- 
eter estimates (treatment effects, standard errors, and variance 
components) by dividing the estimated bias by its true value. 
We considered the cutoff proposed by Hoogland and 
Boomsma (1998) to evaluate the values of relative bias. 
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Values less than .05 are considered as acceptable for the fixed 
effect estimates and the values less than .10 are acceptable for 
the variance component estimates. 

Comparing the results of these analyses in terms of bias, 
MSE, standard errors, and confidence interval coverage pro- 
portions can clarify the added value of performing a multivar- 
iate multilevel meta-analysis over separate univariate multi- 
level meta-analyses for each kind of regression coefficients. 
We did not expect major differences between these methods in 
terms of bias or confidence interval coverage proportions, but 
we did anticipate better precision (i.e., smaller MSE’s and 
standard errors) when performing a multivariate meta- 
analysis as was found in Van den Noortgate et al. (2013). 

Finally, we conducted ANOVA’s to investigate which sim- 
ulation factors affect the bias, MSE, the standard error, and the 
CP. In these analyses, main effect and first order interaction 
effects were included. In our discussion, we only focused on 
effects with eta squared (77°) higher than .01. We chose this 
cutoff because these effects are relatively large compared to 
effects of other design factors, although based on the rule of 
thumb suggested by Cohen (1988) 17° equal to .01 is consid- 
ered as a small effect. 


Results 


Bias and MSE of the overall treatment effect 
estimates 


We looked at the absolute bias of the overall treatment effects 
estimates (i.e., immediate treatment effect and treatment effect 
on time trend). According to an ANOVA, none of the design 
factors had a substantial impact on the bias of overall treat- 
ment effects estimates. We looked at the relative bias of the 
immediate effect estimates and it was close to zero across all 
models. The highest relative bias was found for Model 1 when 
the overall population effect is 2, resulting in a relative bias of 
0.0010. The relative bias for Model 2 and Model 3 was 0.0008 
and 0.0007 respectively. 

We also looked at the MSE of the average treatment effect 
estimates, as an indicator of precision providing important 
information about bias and variance of the estimates. The 
smaller the values of MSE, the better the estimates. The 
ANOVA revealed that the number of studies at the highest 
level (1? = .050) and the between-study variance (1? = .046), 
were the design factors with substantial impact on the MSE of 
the estimated treatment effect and to a smaller extent the in- 
teraction of these factors (17° = .013). An overview of how the 
estimated MSE values is associated with the relevant simula- 
tion factors is shown in Table 2. There were no substantial 
differences between the underlying models in terms of MSE 
estimates. As we can see in Table 2, the largest MSE were 
obtained in conditions with large between-study variance and 
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Table 2. Mean squared error (MSE) of the estimated immediate effect 
(Y200 ) 
on K=10 K= 20 K= 40 
Model 1 0.5 0.196 0.097 0.049 
8 0.953 0.478 0.239 
Models 2 and 3 0.5 0.191 0.095 0.048 
8 0.948 0.476 0.238 


small number of studies. As the number of studies increased 
and the between-study variance decreased, the MSE estimates 
dropped dramatically. As the number of studies increased 
from 10 to 40, the values of MSE became four times smaller. 
The MSE estimates range from 0.016 up to 1.232 with the 
largest value when ten measurement occasions, three cases, 


ten studies, a large value of between-case variance (i.e., a, 


= 8 ), and a large value of between-study variance (i.e., ae 
= 8 ) are included. 


K = number of studies; 02 


= between-study variance; 
Model | = univariate three-level meta-analysis; Model 2 = 
multivariate three-level meta-analysis ignoring covariance; 
Model 3 = multivariate three-level meta-analysis modeling 


covariance. 


Relative standard error bias and confidence interval 
coverage proportion 


An ANOVA indicated that the relative standard error bias was 
mostly affected by the underlying model (77 = .296) and the 
number of studies meta-analyzed (1° = .259). Negative values 
of the relative bias indicate that the standard error estimates 
were slightly smaller than expected. The results in Table 3 
reveal that in none of the conditions, the relative bias exceeds 
the .10 or -.10 cut-off score. The standard errors are more 
underestimated when the sandwich estimator procedure is ap- 
plied, particularly when the number of studies is small. The 
results also indicate as the number of studies increases to 40, 


Table 3 Relative difference between the median of the standard error 
estimates and the standard deviation of the estimated effects 


K=10 K=20 K=40 
Model 1 — 0.024 — 0.009 — 0.003 
Model 2 — 0.022 — 0.008 — 0.002 
Model 2 with RVE — 0.063 — 0.023 — 0.008 
Model 3 — 0.022 — 0.009 — 0.003 
Model 3 with RVE — 0.063 — 0.023 — 0.008 


Notes. K = number of studies; Model | = univariate three-level meta- 
analysis; Model 2 = multivariate three-level meta-analysis ignoring co- 
variance; Model 3 = multivariate three-level meta-analysis modeling 
covariance 


the relative bias is close to zero and negligible across all anal- 
ysis models no matter how many cases and how many mea- 
surements are included. 

We also looked at the confidence interval coverage propor- 
tions to evaluate the estimated standard errors. The estimated 
coverage proportions range from .93 to .96. As can be seen in 
Table 4, there are slight differences in coverage proportion of 
the 95% confidence intervals between different approaches. In 
general, the CPs are close to the nominal values in most of the 
conditions across different models, except in models applying 
sandwich estimator with a small number of studies (K = 10), 
where the CP is slightly below the nominal value due to higher 
underestimation of standard errors in this condition. 


Variance components estimates 


We evaluated the quality of the variance components esti- 
mates (i.e., the between-study and between-case variance) 
for both the immediate treatment effect and the treatment ef- 
fect on time trend. We only report the results associated to the 
immediate treatment effect because similar results were ob- 
tained. The ANOVA indicated that none of the factors had a 
substantial effect on the bias of the between-study variance 
estimates. However, there was notable variation across condi- 
tions. The median relative bias across all design conditions 
was — 0.059, — 0.052, and — 0.055 for Models 1, 2, and 3, 
respectively. 

Table 5 gives an overview of the results across conditions. 
A negative relative bias in the estimated between-study vari- 
ance was found in all conditions. From this table, we can 
deduce that the relative bias ranges from — 0.017 to — 0.388 
when univariate multilevel meta-analysis applied. This range 
equals to — 0.011 to — 0.390 for the approach that ignored the 
covariances at the case and study level and — 0.017 to — 0.355 
for the approach that modeled these covariances. These results 
reveal that the estimates are unbiased as the number of studies 
increase to 40 and the number of cases is either 4 or 7. 
Additionally, the estimates are almost unbiased when the 


Table 4 Coverage proportion for the estimated immediate effect (Y99 ) 


K=10 K=20 K=40 
Model | 0.9538 0.9523 0.9506 
Model 2 0.9532 0.9524 0.9515 
Model 2 with RVE 0.9399 0.9448 0.9471 
Model 3 0.9531 0.9519 0.9508 
Model 3 with RVE 0.9397 0.9447 0.9471 


Notes. K = number of studies; Model 1 = univariate three-level meta- 
analysis; Model 2 = multivariate three-level meta-analysis ignoring co- 
variance; Model 3 = multivariate three-level meta-analysis modeling 
covariance 


Coverage proportions between .94 and .96 are in boldface 
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Table5 Median of relative deviation of the between-study variance estimates of Y299 compared to the true value 


K=10 
z o J=3 J=4 J=7 

Modell 0.5 0.5 -0.188  -0.163  -0.114 
8 —0.087  -0.082 -0.082 

8 0.5 -0.388  -0.360 -0.233 

8 -0.114  -0.100  —0.092 

Model2 0.5 0.5 -0.137 -0.130 — — 0.099 
8 —0.079  -0.080  —0.079 

8 0.5 -0.390  -0.378 0.225 

8 -0.110  —0.099  — — 0.090 

Model3 0.5 0.5 -0.138  -0.137  —0.107 
8 -0.079 0.082  -0.081 

8 0.5 -0.321 -0355 -0.224 

8 -0.105  -0,095 —0,092 


2 


ug 


Notes. J = number of cases; K = number of studies; o 


K=20 K=40 

J=3 J=4 J=7 J=3 J=4 J=7 
— 0.093 — 0.069 — 0.055 — 0.042 — 0.039 —0.021 
— 0.048 — 0.038 — 0.039 —0.019 —0.019 —0.017 
— 0.174 — 0.174 — 0.124 — 0.120 — 0.079 — 0.058 
— 0.064 — 0.050 — 0.045 — 0.026 — 0.024 —0.019 
— 0.064 —0.051 —0.041 — 0.027 — 0.027 — 0.011 
— 0.044 — 0.034 — 0.036 —0.013 —0.014 — 0.014 
— 0.172 — 0.174 —0.121 — 0.120 — 0.083 — 0.060 
— 0.062 — 0.048 — 0.042 — 0.022 —0.021 — 0.017 
— 0.073 — 0.062 — 0.050 — 0.039 — 0.039 —0.021 
— 0.047 — 0.038 — 0.039 —0.019 —0.019 —0.017 
— 0.154 — 0.167 — 0.122 — 0.123 — 0.085 — 0.059 
— 0.064 — 0.050 — 0.045 — 0.025 — 0.025 —0.019 


= between-case variance; 07, = between-study variance; Model 1 = univariate three-level meta- 


analysis; Model 2 = multivariate three-level meta-analysis ignoring covariance; Model 3 = multivariate three-level meta-analysis modeling covariance 


The values in boldface are larger than the .10 cutoff proposed by Hoogland and Boomsma (1998) 


between-study variance is larger (1.e., a, = 8) independent of 
the number of cases, the number of studies, and the value of 
between-case variance. The smallest values for the relative 
bias were found in conditions with a larger number of studies 
(K = 40), independent of other design factors. Across all 
models, the smaller relative bias was obtained when the 
between-study variance is larger. 

The ANOVA showed that the MSE of between-study var- 
lance estimates are mainly affected by the number of studies 
(7° = .021) and the between-study variance (7° = .033), and to 
a minor extent by the interaction of these two factors (7° = 
.011). Table 6 displays the MSE estimates in terms of these 
effective design factors. The MSE estimates are particularly 
large in conditions with a larger between-study variance (1.e., 
g,,, = 8). Closer inspection of the results indicates that the 
MSE decreases significantly as the number of studies K in- 
creases and the between-case and between-study variances 
decrease. The largest MSE was obtained in K = 10, J/=4, J 
= 10, on, = 8, o = 8 condition. We also found that the 
smallest MSE occurs in conditions with a large number of 
studies, a large number of cases, more measurement occa- 
sions, and a small between-case and between-study variabili- 
ty. The results also reveal very similar MSE estimates across 
different approaches. This means that there is no clear advan- 
tage of modeling the covariances at higher levels, in terms of 
the precision of the estimated between-study variance. 

We also looked at the relative bias in the estimate of the 
between-case variance of the immediate treatment effect. The 
preliminary ANOVA showed that the number of measurement 
occasions (0? = .098), the between-case variance (1? = .052), 
and the interaction of between-study variance and the number 
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of measurement occasions (7 = .040) are the most influential 
design factors explaining variability in between-case variance 
estimates. As shown in Table 7, the relative bias exceeds the 
cut-off of .10 when the population between-case variance is 
smaller (i.e., 02, = 0.5) regardless the values of number of 
measurement occasions and between-study variance. 
Specifically, the relative bias is highly larger in condition with 
smaller number of measurement occasions (J = 10). 
Furthermore, the relative bias estimates are smaller in 


Table 6 Mean squared error (MSE) of the estimated between-study 
variance 


- o, K=10 K=20 K= 
40 

Model 1 0.5 0.5 0.220 0.114 0.059 
8 17.545 8.185 4.013 

8 0.5 1.635 0.791 0.414 

8 25.844 12.308 6.091 

Model 2 0.5 0.5 0.209 0.107 0.055 
8 17.462 8.137 4.004 

8 0.5 1.598 0.775 0.406 

8 25.794 12.293 6.084 

Model 3 0.5 0.5 0.215 0.106 0.054 
8 17.626 8.099 3.972 

8 0.5 1.690 0.777 0.403 

8 26.029 12.236 6.037 


Notes. K = number of studies; ny = between-case variance; , = 
between-study variance; Model 1 = univariate three-level meta-analysis; 
Model 2 = multivariate three-level meta-analysis ignoring covariance; 


Model 3 = multivariate three-level meta-analysis modeling covariance 
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Table 7 Median of relative deviation of the between-case variance es- 
timates of Y299 


z a T=10 1=20 1=40 
Model 1 0.5 0.5 0.271 0.191 0.094 
8 1.928 0.692 0.308 
8 0.5 0.083 0.017 -0.005 
8 0.213 0.065 0.022 
Model 2 0.5 0.5 2.083 0.532 0.177 
8 3.760 1.045 0.393 
8 0.5 0.215 0.046 0.002 
8 0.347 0.095 0.029 
Model 3 0.5 0.5 2.102 0.508 0.170 
8 3.762 1.013 0.385 
8 0.5 0.201 0.037 —0.000 
8 0.331 0.085 0.026 


2 


‘, = between-case vari- 


Notes. I = number of measurement occasions; 0 
ance; a. = between-study variance; Model 1 = univariate three-level 
meta-analysis; Model 2 = multivariate three-level meta-analysis ignoring 
covariance; Model 3 = multivariate three-level meta-analysis modeling 


covariance 


The values in boldface are larger than the .10 cutoff proposed by 
Hoogland and Boomsma (1998) 


conditions with larger between-case variance and a number of 
measurement occasions larger than 10 (i.e., J = 20 or 40). In 
general, the univariate approach led to better estimation of the 
between-case variance compared to other approaches. The 
approaches of ignoring and modeling the covariances in sec- 
ond and third levels obtained very similar results. 

According to the ANOVA, there were no design factors 
that substantially affected the MSE for between-case variance 
estimates. The highest 17° was observed for the number of 
measurement occasions (7° = .004). As Table 8 shows, the 
MSE is noticeably large when the number of measurement 
occasions is 10. The MSE’s highly drop as the number of 
measurement occasions increases. As was the case for the 
between-study variance, the MSE for the between-case vari- 
ance was smallest when the number of observations, cases and 
studies was largest and the variances at the higher levels 
smaller. 


Bias and MSE of between-study covariance 
and between-case covariance 


The overall absolute bias of the estimated between-case co- 
variance (1.¢., 0213) equals 0.024 and ranges from — 0.008 to 
0.090 depending on the condition. We did not look at relative 
bias for covariances because of having value of 0 for covari- 
ances and it was not possible to calculate relative bias for some 
conditions. According to the ANOVA, the number of mea- 
surement occasions (f° = .013) is related to the bias (with a 
higher / corresponding to less bias). The MSE values ranged 


from zero to 0.679. The largest MSE was observed in K = 10, 
J=4,/= 10, Cn, ='8, o., = 8 condition. Similar to the bias, 
the number of measurement occasions (17° = .006) appeared to 
be the factor most influencing variability of MSE estimates. 
As can be seen in Fig. 1, the bias and MSE is greatly reduced 
as the number of measurement occasions increases from 10 to 
20. 

The mean bias of the between-study covariance was equal 
to — 0.003 with the range between — 0.023 and 0.007 
representing unbiased estimates. For the estimated between- 
study covariance, no substantial effects of the design factors 
on the bias were found. The MSE values ranged from 0.006 to 
0.837 with the largest value in K = 10, /=4, J= 10, a. = 5: 
o;, = 8 condition. The ANOVA showed that the number of 
studies (1° = .032), the number of measurement occasions (0? 
= .016), and between-study variance (17° = .016) have signif: 
icant effects on MSE estimates. Table 9 indicates that the 
smallest MSE values were observed in condition with the 
number of measurement occasions larger than 10 (i.e., 7 = 20 
or 40) and the smaller between-study variance (i.e., a = 0.5). 

In sum, the results indicate that the between-case and 
between-study covariance are well estimated in terms of ac- 
curacy and precision. Also estimating the between-case and 
especially the between-study variance, results in small and 
even negligible bias in almost all conditions. Our results also 
indicate that the approach modeling the covariances does not 
significantly outperform the approach ignoring the covari- 
ances at the case and study level. 


Discussion 


In the current simulation study, we studied the meta-analysis 
of multiple regression coefficients from SCEDs. We simulated 
SCED raw data with the covariance between the immediate 
treatment effect and the treatment effect on the time trend at 
the case and study level. Afterwards, we used five main ap- 
proaches to handle the existing dependence between the main 
treatment effects (i.e., the immediate treatment effect and the 
treatment effect on the time trend). 


Table 8 Mean squared error (MSE) of the estimated between-case 
variance 


T=10 [=20 1=40 
Model 1 10.850 2.241 1.527 
Model 2 13.947 2.437 1.544 
Model 3 13.586 2.344 1.531 


Notes. = number of measurement occasions; Model | = univariate three- 
level meta-analysis; Model 2 = multivariate three-level meta-analysis 
ignoring covariance; Model 3 = multivariate three-level meta-analysis 
modeling covariance 
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Fig. 1 Bias and MSE of the estimated between-case covariance esti- 
mates; / = number of measurement occasions 


In the first approach, we ignored the dependence among 
the treatment effects by running a separate univariate three- 
level meta-analysis for each treatment effect. In the second 
approach, we used a multivariate model but ignored the 
existing simulated covariance at the case and study level. In 
the third approach, we used the same model but used a sand- 
wich estimator for the standard errors of the effect estimates. 
The fourth approach used again a multivariate model but in- 
cluded the covariance between the treatment effects at higher 
levels. In the fifth approach, a sandwich estimator was applied 
for the standard errors for this multivariate model. The ques- 
tion was whether handling the dependency among multiple 
regression coefficients for each case within SCED studies in 
a three-level meta-analysis with these approaches affects the 
estimation of statistical properties of fixed effects and variance 
components. 

Based on the reported results, we did not find substantial 
differences between the underlying models in terms of the bias 
and MSE of the treatment effect estimates. However, the uni- 
variate model (i.e., Model 1) led to slightly larger estimates of 
the overall treatment effects compared to other models. None 
of the simulation design factors had a substantial impact on the 
treatment effects estimates. In contrast, the number of studies, 
the between-study and the between-case variance were found 


Table 9 Mean squared error (MSE) of the estimated between-study 
covariance 


I o, K=10 K=20 K=40 
10 0.5 0.281 0.171 0.130 
8 0.443 0.301 0.227 
20 0.5 0.092 0.081 0.077 
8 0.173 0.158 0.151 
40 0.5 0.081 0.076 0.073 
8 0.161 0.153 0.148 


Note. = number of measurement occasions; K = number of studies; Oo, = 
between-study variance 
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to have significant effects on MSE. The highest MSE values 
were obtained in the conditions with a larger between-study 
variance and a small number of studies. In terms of the stan- 
dard errors of the treatment effects, the underlying models 
along with the number of studies and the interaction of num- 
ber of measurement occasions, number of cases, and number 
of studies played a significant role in the variability of the 
estimates. The results indicate a greater underestimation when 
the sandwich estimator was applied to results from both ap- 
proaches of modeling and when ignoring the covariances at 
higher levels. This in turn led to smaller coverage proportions 
for results from the models paired with the sandwich estima- 
tor, which was slightly smaller than the nominal coverage 
proportion values. In general, the overall treatment effects 
are estimated precisely and accurately when the covariance 
exists at the case and study level, independent of the underly- 
ing model fitted. The results indicate that the estimates of the 
treatment effects are robust even when the random effects in 
the model are misspecified through ignoring the existing 
covariances. 

The between-study and between-case variances are not es- 
timated precisely and accurately in certain conditions. We 
found underestimated between-study variances across all 
combinations of conditions. However, unbiased estimates 
were obtained for the between-study variance in most condi- 
tions. Biased between-study variance estimates were observed 
particularly in scenarios with 10 or 20 studies and large 
between-case variances. On the other hand, more bias was 
observed in the between-case variance estimates. Conditions 
with a larger between-case variance and 20 or 40 measure- 
ment occasions resulted in unbiased between-case variance 
estimates. Biased estimates were obtained in other conditions. 
In general, modeling the covariances at higher levels led to 
slightly less biased estimations of the between-case and 
between-study variances in conditions with a larger 
between-case variance. The reported results indicate that the 
approach that models the existing covariances does not out- 
perform the approach that ignores these covariances. Also, in 
this simulation, we obtained biased and less precise estimates 
for the between-case variance compared to the estimates for 
the between-study variance. 

These results are not completely in line with the findings in 
the study of Moeyaert et al. (2016), who found more bias for 
the between-study variance than the between-case variance. 
An explanation is that we simulated the covariance between 
the immediate treatment effect and the treatment effect on the 
time trend at the higher levels, whereas Moeyaert et al. (2016) 
generated the covariance between the baseline level and the 
treatment effect not between two dependent treatment effect 
sizes. Our results show that the degree of bias does not differ 
substantially across the models (1.e., the misspecified model 
that ignores the covariances at higher levels and the correctly 
specified model), and unbiased estimates are obtained for the 
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between-case and the between-study covariance in different 
conditions. 


Limitations and recommendations 


The current study suffers from some limitations. The general- 
izability of the findings in this study is one of the major con- 
cerns, although we generated the data considering certain re- 
alistic conditions that seem representative for the single-case 
experimental studies. We limited the simulated three-level da- 
ta to data from MBDs. We also did not include a possible 
autocorrelation between the residuals at the first level although 
the existence of autocorrelation might be expected when a 
case is repeatedly measured. Furthermore, we only included 
the correlation between the regression coefficients associated 
with the treatment effects and set the correlation between other 
coefficients to zero to explore the performance of different 
models in simple scenarios. We only generated conditions 
with the covariances at the case and study level together. We 
did not simulate conditions with the covariances at only one of 
these levels. 

Moreover, in this study we evaluated the performance of 
different models in terms of handling dependence among the 
treatment effects at higher levels applying the three-level me- 
ta-analytic approach. This is only one possible source of de- 
pendence among effect sizes. Dependence among effect sizes 
might occur due to several other reasons. Multiple effect sizes 
might be obtained by comparison of different treatments, by 
using multiple measures for the same construct or by using 
multiple outcome variables (Cheung, 2014; Van den 
Noortgate et al., 2013). Dependent effect sizes are less infor- 
mative compared to independent effect sizes because they are 
calculated from common or related data and therefore give 
partly the same information. Thus, including dependent effect 
sizes in the meta-analysis as independent effects leads to using 
the same information multiple times that in turn, would over- 
estimate confidence in the results of the meta-analysis. By 
treating dependent effect sizes this way, an underestimation 
of the standard errors is likely, resulting in too small confi- 
dence interval (Van den Noortgate et al., 2013). It might be 
interesting to investigate other sources of dependence among 
effect sizes from SCED studies in future studies. 

This simulation further generated data by sampling resid- 
uals from normally distributed and homogeneous population 
distributions, and assuming that time trends are linear. Finally, 
we meta-analyzed only treatment effect sizes calculated as 
regression coefficients. It might be interesting if in future re- 
search other parametric and non-parametric effect sizes 
(Manolov & Moeyaert, 2017) are synthesized applying a mul- 
tilevel meta-analytic approach. Some non-overlap indices 
such as percentage of non-overlapping data (PND; Parker, 
Vannest, & Davis, 2011; Scruggs, Mastropieri, & Casto, 
1987), percentage of all non-overlapping data (PAND; 


Parker, Hagan-Burke, & Vannest, 2007; Parker et al., 2011), 
and percentage of zero data (PZD; Campbell, 2004) or a stan- 
dardized mean difference (SMD; Gingerich, 1984; Maggin 
et al., 2011) can be calculated instead of regression coeffi- 
cients as the treatment effect sizes. 

In summary, the findings in the current study indicate no 
substantial differences between the results of applying differ- 
ent three-level meta-analytic models, confirming the robust- 
ness of the three-level approach in handling the dependence 
among the treatment effects at higher levels even if the covari- 
ance structure is misspecified. Although biased variance com- 
ponents estimates were obtained in some conditions, the 
misspecification of the covariance structure did not lead to 
worse estimates. We recommend SCED researchers use a 
multivariate multilevel approach to model the existing depen- 
dence among multiple regression coefficients (i.e., treatment 
effect sizes). Using this approach gives the opportunity to 
model the covariance at higher levels when combining the 
results from the primary SCEDs. 
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